
rm(list = ls())

library(magrittr)
library(ggplot2)

### Read in data
stock   <- read.csv(paste0("C:/Users/tbrai/Dropbox/DBT/Data/ePOS RCT Jharkhand/New Data Stream/Adherence/ForPlots/Data/adherence_disbursementstock_compare2.csv"))

### Prep data 
stock              <- plm::pdata.frame(stock, index=c("dealer_id_num", "month"))
stock$agg_stock_l1 <- lag(stock$agg_stock, k = 1)
stock$lag_carry    <- lag(stock$agg_carryover, k = 1)

### Gen vars for calculated stock 
stock$stock_calculated <- stock$agg_stock_l1 + stock$dis_qty - stock$agg_offtake

### Gen vars for calculated disbursement
stock$dis_calculated <- stock$ent_qty + stock$agg_carryover - stock$agg_stock_l1
stock$dis_calculated[stock$dis_calculated < 0] <- 0

### Calculated stock indexed at month 1
stock <- dplyr::group_by(stock, dealer_id_num)
stock <- dplyr::mutate(stock, stock_calc_t1 = dplyr::first(agg_stock))

stock <- dplyr::mutate(stock, cum_dis = ifelse(is.na(dis_qty), 0, dis_qty))
stock <- dplyr::mutate(stock, cum_dis = cumsum(cum_dis))
stock <- dplyr::mutate(stock, cum_off = ifelse(month == 0, 0, agg_offtake))
stock <- dplyr::mutate(stock, cum_off = cumsum(cum_off))
stock <- dplyr::mutate(stock, cum_stock = stock_calc_t1 + cum_dis - cum_off)

### Panel 1: Measures of stock
stocks_plot <- stock %>% 
  dplyr::filter(month != 0) %>%
  dplyr::mutate(month = as.numeric(month) - 1) %>% 
  dplyr::group_by(month) %>%
  dplyr::summarize(mean_stock   = mean(agg_stock, na.rm = TRUE), 
                   mean_stock2  = mean(cum_stock, na.rm = TRUE), 
                   mean_ent     = mean(ent_qty, na.rm = TRUE)) %>% 
  tidyr::pivot_longer(cols = c(mean_stock, mean_stock2, mean_ent)) %>% 
  ggplot(., mapping = aes(x = month, y = value, color = name)) + 
  geom_line(size = 1.5, alpha = 0.8) +
  scale_x_continuous(breaks = rep(1:11), 
                 labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov")) +
  scale_y_continuous(breaks = seq(0, 25000, 2500)) + 
  scale_color_manual(breaks = c("mean_stock", "mean_stock2", "mean_ent"), 
                     values = c("black", "#4669bd", "grey70"), 
                     labels = c("Recorded stock", "Counterfactual stock", "Entitlement")) +
  labs(x = "", y = "", color = "", title = "Panel B: Measures of stock") + 
  theme_bw() + 
  theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text = element_text(size = 20), 
        axis.title = element_text(size = 20), 
        legend.text = element_text(size = 20),
        title = element_text(size = 16))
stocks_plot
ggsave(plot = stocks_plot, filename = paste0(here::here(), "/exhibits/all_possible_stock_measures.png"))

### Panel 2: Calculated disbursement using cumulative stock
stock                  <- dplyr::group_by(stock, dealer_id_num)
stock                  <- dplyr::mutate(stock, cum_stock_l1 = dplyr::lag(cum_stock, n = 1))
stock                  <- dplyr::mutate(stock, cum_disbursement = ifelse(as.numeric(month) <= 7, ent_qty, ent_qty + lag_carry - cum_stock_l1))
stock                  <- dplyr::mutate(stock, cum_disbursement = ifelse(cum_disbursement < 0, 0, cum_disbursement))

dis_plot <- stock %>% 
  dplyr::filter(month != 0) %>% 
  dplyr::mutate(month = as.numeric(month) - 1) %>% 
  dplyr::group_by(month) %>%
  dplyr::summarise(mean_cum_dis  = mean(cum_disbursement, na.rm = TRUE),
                   mean_dis_qty  = mean(dis_qty, na.rm = TRUE),
                   mean_ent_qty  = mean(ent_qty, na.rm = TRUE)) %>%
  tidyr::pivot_longer(cols = c(mean_cum_dis, mean_dis_qty, mean_ent_qty)) %>% 
  ggplot(., mapping = aes(x = month, y = value, color = name)) + 
  geom_line(size = 1.5, alpha = 0.8) +
  scale_x_continuous(breaks = rep(1:11), 
                     labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov")) +
  scale_color_manual(breaks = c("mean_dis_qty", "mean_cum_dis", "mean_ent_qty"), 
                     values = c("black", "#4669bd", "grey70"), 
                     labels = c("Recorded disbursement", "Counterfactual disbursement", "Entitlement")) +
  scale_y_continuous(breaks = seq(0, 7000, 1000),
                     limits = c(0, 7000)) +
  labs(x = "Month", y = "", color = "", title = "Panel B: Measures of disbursement") + 
  theme_bw() + 
  theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text = element_text(size = 20), 
        axis.title = element_text(size = 20), 
        legend.text = element_text(size = 20),
        title = element_text(size = 16))
dis_plot
ggsave(plot = dis_plot, filename = paste0(here::here(), "/exhibits/all_possible_disbursement_measures.png"))

combined <- gridExtra::grid.arrange(gridExtra::arrangeGrob(stocks_plot, dis_plot, ncol = 1, 
                                                           left = grid::textGrob("Value (kg)", rot = 90,
                                                                                 gp = grid::gpar(fontsize = 14))))
ggsave(plot = combined, filename = paste0(here::here(), "/replication/exhibits/stock_disbursement_comparison.pdf"), width = 15, height = 12)

### Summary statistics related to the above graph
mean_stock_by_month <- stock %>% 
  dplyr::filter(month != 0) %>%
  dplyr::mutate(month = as.numeric(month) - 1) %>% 
  dplyr::group_by(month) %>%
  dplyr::summarize(mean_stock   = mean(agg_stock, na.rm = TRUE), 
                   mean_stock2  = mean(cum_stock, na.rm = TRUE), 
                   mean_ent     = mean(ent_qty, na.rm = TRUE))
  
rec_diff <- abs(mean_stock_by_month$mean_stock[mean_stock_by_month$month == 9] - mean_stock_by_month$mean_stock[mean_stock_by_month$month == 6])

calc_diff <- abs(mean_stock_by_month$mean_stock2[mean_stock_by_month$month == 9] - mean_stock_by_month$mean_stock2[mean_stock_by_month$month == 6])

t.test(stock$agg_stock, stock$agg_stock)

t.test(abs(stock$agg_stock[stock$month == 9] - stock$agg_stock[stock$month == 6]),
       abs(stock$cum_stock[stock$month == 9] - stock$cum_stock[stock$month == 6]))

t.test(stock$agg_stock - stock$cum_stock, mu = 0)

t.test(abs(stock$agg_stock[stock$month == 9] - stock$agg_stock[stock$month == 6]) -
       abs(stock$cum_stock[stock$month == 9] - stock$cum_stock[stock$month == 6]), mu = 0)
